{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 4 : Calculus and power series\n",
    "**References:**\n",
    "* [Introduction to Symbolic Computation](http://homepages.math.uic.edu/~jan/mcs320/mcs320notes/index.html), Lectures \n",
    "[18](http://homepages.math.uic.edu/~jan/mcs320/mcs320notes/lec18.html), [19](http://homepages.math.uic.edu/~jan/mcs320/mcs320notes/lec19.html), [21](http://homepages.math.uic.edu/~jan/mcs320/mcs320notes/lec21.html), [22](http://homepages.math.uic.edu/~jan/mcs320/mcs320notes/lec22.html), [23](http://homepages.math.uic.edu/~jan/mcs320/mcs320notes/lec23.html)\n",
    "* [Calculus in SageMath for high-school students](https://www.sagemath.org/calctut/index.html)\n",
    "* [Calculus in SageMath for first-year university students](https://doc.sagemath.org/html/en/prep/Calculus.html)\n",
    "* [Sage Quickstart for Multivariable Calculus](https://doc.sagemath.org/html/en/prep/Quickstarts/Multivariable-Calculus.html)\n",
    "* [Vector calculus in SageMath](https://doc.sagemath.org/html/en/thematic_tutorials/vector_calculus.html)\n",
    "* [Power Series rings in SageMath](https://doc.sagemath.org/html/en/reference/power_series/sage/rings/power_series_ring.html)\n",
    "\n",
    "**Summary:**<br>\n",
    "In this lecture, we start by recalling how new **functions** can be defined in Python. Then we show how to **plot** functions in SageMath and discuss concepts from calculus (or analysis) like **limits**, **derivatives** and **integrals**. We finish by discussing **power series**, which are a generalization of the Taylor series from analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python warmup: functions\n",
    "We have already seen many functions implemented in SageMath (such as ``log``, ``span`` or ``solve``). But it is also possible (and often very useful) to write functions ourselves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def parabola(x):\n",
    "    x_squared = x^2\n",
    "    return x_squared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This creates a new function ``parabola``, which we can immediately call on some value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parabola(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The general code pattern for defining a function ``foo`` with $n$ inputs (or **arguments**) looks as follows:\n",
    "```\n",
    "def foo(arg1, arg2, ..., argn):\n",
    "    some_code()\n",
    "    more_code()\n",
    "    return result\n",
    "```\n",
    "As was the case with ``if`` and ``for``, we have a block of indented code (the **body** of the function) in which we can do what we want (e.g. create new variables, call other functions, use ``if`` and ``for``). Most functions ``foo`` should in the end give back some result, which we do using the ``return`` keyword.\n",
    "\n",
    "Let's see a function with more than one argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bigger_than(x,y):\n",
    "    if x > y:\n",
    "        return True\n",
    "    else:\n",
    "        return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bigger_than(5,6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "Define the following functions and test them in some examples.\n",
    "* the function ``sign`` returning the [sign](https://en.wikipedia.org/wiki/Sign_function) of a real number\n",
    "* the function ``conjugate``, taking $n \\times n$-matrices $A,S$ and returning $S A S^{-1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution** (uncomment to see)\n",
    "<!---\n",
    "* the function ``sign``\n",
    "\n",
    "```\n",
    "def sign(x):\n",
    "    if x<0:\n",
    "        return -1\n",
    "    if x==0:\n",
    "        return 0\n",
    "    if x>0:\n",
    "        return 1\n",
    "sign(-34)\n",
    "> -1\n",
    "sign(0)\n",
    "> 0\n",
    "sign(13.4)\n",
    "> 1\n",
    "```\n",
    "* the function ``conjugate``\n",
    "\n",
    "```\n",
    "def conjugate(A, S):\n",
    "    return S * A * S.inverse()\n",
    "A = matrix([[1,2],[3,4]])\n",
    "S = matrix([[0,1],[1,0]])\n",
    "conjugate(A,S)\n",
    "> \n",
    "[4 3]\n",
    "[2 1]\n",
    "```\n",
    "--->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also have functions with **optional arguments**. These are arguments which *can* be given, but don't *have to* be given when calling the function. However, in exchange you need to provide a *default value* for each of them. Let's generalize the function ``parabola`` from above to implement the function $x \\mapsto a \\cdot x^2 + c$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def parabola(x, a=1, c=0):\n",
    "    x_squared = x^2\n",
    "    return a*x_squared + c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parabola(2, 5, 1) # returns 5 * 2^2 +1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parabola(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parabola(3,-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we see above, if a function gets a number of arguments which is smaller than the total number of fixed and optional arguments, it simply fills the optional arguments in the order that was specified in the function definition. \n",
    "\n",
    "If we want to explicitly specify the argument ``c`` above, but leave ``a`` to take the default value, we can do so by explicitly specifying the name of the argument that we want:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parabola(3, c=-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The general pattern looks like this:\n",
    "```\n",
    "def foo(arg1, arg2, ..., argn, opt_arg1 = default1, ..., opt_argm = defaultm):\n",
    "    some_code()\n",
    "    more_code()\n",
    "    return result\n",
    "```\n",
    "Here, it is **important** that the non-optional arguments come before the optional ones."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "Implement a function ``fibonacci`` computing the Fibonacci-Numbers $F_n$ defined by\n",
    "$$\n",
    "F_0 = 0, \\quad F_1 = 1,  \\quad F_{n} = F_{n-1} + F_{n-2} \\text{ for }n \\geq 2.\n",
    "$$\n",
    "Add the option to specify $F_0, F_1$ by hand."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution** (uncomment to see)<br>\n",
    "<!---\n",
    "While it is of course possible to implement this using ``for``-loops and lists or dictionaries, we'll do it with a *recursive* function instead.\n",
    "```\n",
    "def fibonacci(n, F0=0, F1=1):\n",
    "    if n == 0:\n",
    "        return F0\n",
    "    if n == 1:\n",
    "        return F1\n",
    "    return fibonacci(n-1,F0,F1) + fibonacci(n-2,F0,F1)\n",
    "```\n",
    "Then we can play around a bit:\n",
    "```\n",
    "[fibonacci(i) for i in range(10)]\n",
    "> [0, 1, 1, 2, 3, 5, 8, 13, 21, 34]\n",
    "[fibonacci(i,0,2) for i in range(10)]\n",
    "> [0, 2, 2, 4, 6, 10, 16, 26, 42, 68]\n",
    "[fibonacci(i,-1,1) for i in range(10)]\n",
    "> [-1, 1, 0, 1, 1, 2, 3, 5, 8, 13]\n",
    "```\n",
    "--->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculus\n",
    "SageMath can do many computations related to the infinitesimal behaviour of functions, such as limits, derivatives and integrals. In each case, the function should be represented by a symbolic expression, as seen in the last lecture. Since it is useful for illustrating the later computations (and in general), we'll start however with some information how to plot functions (a.k.a. producing pretty pictures)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting\n",
    "The most basic task here is to plot a function in one variable on an interval. For instance, we can give the function as a symbolic expression in ``x``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(sin(x),x,-5,5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can program a function as seen in the first part of the lecture:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def g(x):\n",
    "    if x <= 0:\n",
    "        return cos(x)\n",
    "    else:\n",
    "        return floor(x) # floor is the round-down function, sending 3.14 to 3, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(g,x,-5,5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the time it's more useful to use the symbolic expressions, though. \n",
    "\n",
    "A curiosity: if you type ``plot(g(x),x,-5,5)`` it does not do what you want. If you understand what goes wrong there, you already have a pretty good grasp on symbolic expressions!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "* Plot your favorite function on an interval of your choice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Plot your second-favorite function on the interval $(-3,3)$, but in green.<br>\n",
    "*Hint:* Since unfortunately I forgot to say how to change the color, you will have to find it out yourself!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution** (uncomment to see)<br>\n",
    "<!---\n",
    "Since I don't know your favorite functions, I'll have to draw some of mine. The first is the probability density function of the [normal distribution](https://en.wikipedia.org/wiki/Normal_distribution) around $0$ with standard deviation $1$:\n",
    "$$\n",
    "f(x) = \\frac{1}{\\sqrt{2 \\pi}} e^{-\\frac{1}{2} x^2}\n",
    "$$\n",
    "Here is the code:\n",
    "```\n",
    "plot(1/sqrt(2*pi)*exp(-1/2*x^2), x, -6, 6)\n",
    "```\n",
    "<img src=\"https://johannesschmitt.gitlab.io/mat007/normaldist.png\" width=500/>\n",
    "Another nice function is\n",
    "$$\n",
    "g(x) = \\frac{\\cos(42 x)}{x^2+1} \n",
    "$$\n",
    "which we plot in green:\n",
    "\n",
    "```\n",
    "plot(cos(42*x)/(x^2+1), x, -3, 3, color = 'green')\n",
    "```\n",
    "\n",
    "<img src=\"https://johannesschmitt.gitlab.io/mat007/greenfunction.png\" width=500/>\n",
    "--->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can guess, there are lots more optional arguments for the ``plot`` function, which you can again find in the documentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(sign, x, -3, 3, ymin=-2, ymax=2, axes=False, fill='axis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function ``plot`` actually returns a Graphics-object and you can *literally* add multiple such objects to plot multiple functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = plot(sin(x),x,-6,6)\n",
    "g = plot(sin(2*x),x,-6,6,color='red')\n",
    "f+g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Three more variants of plots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('y')\n",
    "implicit_plot(x^2+y^2-1, (-3,3), (-2,2))  # plots the solution set of x^2 + y^2 -1 = 0 inside [-3,3] x [-2,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('t')\n",
    "parametric_plot((3*cos(t), 2*sin(t)), (t,0,2*pi)) # plots the curve t -> (3 cos(t), 2 sin(t)) for t in [0, 2 pi]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plot3d(sin(x) + e^cos(y),(x,-3,3),(y,-3,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "Consider the following system of equations:\n",
    "$$\n",
    "7 \\, x^{2} - 11 \\, x y + 14 \\, y^{2} + 4 \\, x - 15 \\, y - 51 = 0\\,,\\\\\n",
    "-49 \\, x^{2} + 59 \\, x y - 94 \\, x - 43 \\, y + 159 = 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Plot the two solution sets to the two equations in separate colors (in the same coordinate system).\n",
    "* Solve this system of equations (using the methods from last lecture) and verify graphically that you found all solutions in the first part of the exercise.\n",
    "\n",
    "To spare you the typing, here are the left-hand sides of the equations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('x,y')\n",
    "F1 = 7*x^2 - 11*x*y + 14*y^2 + 4*x - 15*y - 51\n",
    "F2 = -49*x^2 + 59*x*y - 94*x - 43*y + 159"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution** (uncomment to see)\n",
    "<!---\n",
    "```\n",
    "A = implicit_plot(F1, (x,-5,6), (y,-5,5))\n",
    "B = implicit_plot(F2, (x,-5,6), (y,-5,5), color='red')\n",
    "A+B\n",
    "```\n",
    "<img src=\"https://johannesschmitt.gitlab.io/mat007/F1F2solution.png\" width=500/>\n",
    "\n",
    "```\n",
    "solve([F1==0, F2==0],[x,y])\n",
    "> [[x == 2, y == 3], [x == 1, y == -1], [x == -1, y == 2], [x == -3, y == 0]]\n",
    "```\n",
    "--->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Limits\n",
    "For a function given by a symbolic expression, we can compute its limit at a given point as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "limit((x^2 - 1)/(x - 1), x=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "Compute the limit\n",
    "$$\n",
    "\\lim_{x \\to 0} \\frac{\\sinh(x)}{x}\n",
    "$$\n",
    "and check with a plot that this is plausible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution** (uncomment to see)\n",
    "<!---\n",
    "```\n",
    "limit(sinh(x)/x, x=0)\n",
    "> 1\n",
    "plot(sinh(x)/x, x, -2, 2)\n",
    "```\n",
    "<img src=\"https://johannesschmitt.gitlab.io/mat007/sinh.png\" width=500/>\n",
    "--->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also deal with divergent functions, one-sided limits and limits as $x$ goes to $\\pm \\infty$: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "limit(1/x, x=0, dir='+')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "limit((x+1)/x, x=Infinity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Derivatives and integrals\n",
    "Of course SageMath can differentiate (symbolic) functions, and for many functions it can also integrate them. Here are three ways to obtain the derivative:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = x * exp(x)\n",
    "print(diff(f,x))\n",
    "print(derivative(f,x))\n",
    "print(f.derivative(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the derivative at a specific point, we can evaluate the symbolic expression there, as seen in the last lecture:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = f.derivative()\n",
    "g(x=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "Write a function ``plot_tangent`` which takes as input\n",
    "* a symbolic expression ``f`` of a function in the variable ``x``\n",
    "* a point ``x0``\n",
    "\n",
    "and outputs a plot of the function ``f`` on the interval $[x_0-3, x_0+3]$ together with its tangent line at ``x0``.<br>\n",
    "*Bonus:* If you are done early, feel free to include optional parameters to change the interval of plotting, the colour of the tangent line etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution** (uncomment to see)<br>\n",
    "<!---\n",
    "Here is some code for the basic function allowing us to plot. Note that I inserted some comments using the symbol ``#``. This is a good idea for more complicated code, both to make it readable to others (and your future self) and for structuring your thoughts.\n",
    "```\n",
    "def plot_tangent(f,x0):\n",
    "    # computing the equation of the tangent line L (as a function in x)\n",
    "    m = f.derivative(x)(x=x0)  # slope of the tangent\n",
    "    L = m*(x-x0) + f(x=x0)\n",
    "    \n",
    "    # plot of f and L\n",
    "    fplot = plot(f,x,x0-3,x0+3)\n",
    "    Lplot = plot(L,x,x0-3,x0+3, color='red')\n",
    "    return fplot + Lplot\n",
    "```\n",
    "Here is an example of an application of our new function:\n",
    "```\n",
    "plot_tangent(sin(x), pi/2+0.1)\n",
    "```\n",
    "<img src=\"https://johannesschmitt.gitlab.io/mat007/tangent.png\" width=500/>\n",
    "--->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course we can also compute partial derivatives for functions depending on multiple parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('a,b')\n",
    "derivative(a^2*b^3,a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is an example of how to do indefinite integration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = integral(sin(x)*cos(x),x); g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The syntax for a definite integral like\n",
    "$$\n",
    "\\int_0^{2 \\pi} g(x) dx\n",
    "$$\n",
    "is as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integral(g, (x,0,2*pi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To apply the things we learned, here is a task from an actual [Analysis II-exam](https://www2.math.ethz.ch/education/bachelor/lectures/fs2016/other/analysis2_mavt_matl/pruefungen/pruefung_winter_16.pdf) at ETH Zurich (see Exercise 3).\n",
    "#### Exercise\n",
    "Let $a, b > 0$. Let the function $f : \\mathbb{R}^2 \\to \\mathbb{R}$ be defined by\n",
    "$$\n",
    "f(x,y) := (a x^2 + b y^2) \\cdot e^{-(x^2 + y^2)}.\n",
    "$$\n",
    "Find the global maxima and minima of $f$ for\n",
    "* $a>b$.\n",
    "* $a=b=1$.\n",
    "\n",
    "*Hint:* The function $f$ slopes down to zero towards infinity.\n",
    "\n",
    "In case you have forgotten how these extreme value problems work, you can of course [look it up on wikipedia](https://en.wikipedia.org/wiki/Maxima_and_minima#Functions_of_more_than_one_variable)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution** (uncomment to see)<br>\n",
    "<!---\n",
    "To find candidates for the local maxima and minima, we use that the partial derivatives of $f$ must be zero at these points:\n",
    "```\n",
    "var('a,b,x,y')\n",
    "f = (a*x^2+b*y^2)*exp(-(x^2+y^2))\n",
    "fx = f.derivative(x)\n",
    "fy = f.derivative(y)\n",
    "solve([fx==0, fy==0], [x,y])\n",
    "> [[x == 0, y == 0], [x == -1, y == 0], [x == 1, y == 0], [x == 0, y == -1], [x == 0, y == 1]]\n",
    "```\n",
    "Then we should substitute these points into the function to see the values of $f$ there. For this, it is more convenient to obtain the solutions as a dictionary:\n",
    "```\n",
    "solution_list = solve([fx==0, fy==0], [x,y], solution_dict=True); solution_list\n",
    "> [{x: 0, y: 0}, {x: -1, y: 0}, {x: 1, y: 0}, {x: 0, y: -1}, {x: 0, y: 1}]\n",
    "for sol in solution_list:\n",
    "    print((sol[x], sol[y]), f.subs(sol))\n",
    ">\n",
    "(0, 0) 0\n",
    "(-1, 0) a*e^(-1)\n",
    "(1, 0) a*e^(-1)\n",
    "(0, -1) b*e^(-1)\n",
    "(0, 1) b*e^(-1)\n",
    "```\n",
    "From its formula, it's clear that $f$ is non-negative, and the unique zero is at $(0,0)$. In fact, the non-negativity can be checked with SageMath:\n",
    "```\n",
    "assume([a>0, b>0])\n",
    "bool(f>=0)\n",
    "> True\n",
    "```\n",
    "This shows that $(0,0)$ is the unique global minimum. Since $f$ goes to zero when $(x,y)$ goes to infinity, the global maxima are precisely the local maxima with the largest values of $f$. Since $a > b$ by assumption, we see that they are exactly at $(\\pm 1, 0)$ with value $a e^{-1}$.\n",
    "\n",
    "Let's repeat the process for $a=b=1$. Here we see that there are solutions which the program above did not find:\n",
    "```\n",
    "f_sp = f.subs({a:1, b:1})\n",
    "f_sp_x = f_sp.derivative(x)\n",
    "f_sp_y = f_sp.derivative(y)\n",
    "sp_solution_list = solve([f_sp_x==0, f_sp_y==0], [x,y], solution_dict=True); sp_solution_list\n",
    "> \n",
    "[{x: r1, y: sqrt(-r1^2 + 1)},\n",
    " {x: r2, y: -sqrt(-r2^2 + 1)},\n",
    " {x: 0, y: 0},\n",
    " {x: -1, y: 0},\n",
    " {x: 1, y: 0}]\n",
    "```\n",
    "If you stare at this a bit, you see that the first two solutions exactly parameterize the upper and lower part of the unit circle using formal parameters ``r1, r2``. If we plot this function, we can see why this is the case:\n",
    "```\n",
    "plot3d(2*f_sp, (x,-2,2), (y,-2,2)) # I multiply the function by 2 for dramatic effect\n",
    "```\n",
    "<img src=\"https://johannesschmitt.gitlab.io/mat007/f_sp_plot.png\" width=800/>\n",
    "\n",
    "We see that the function has a (continuous) maximum along the unit circle. We also see this by plugging in the solutions we found, as before:\n",
    "```\n",
    "for sol in sp_solution_list:\n",
    "    print((sol[x], sol[y]), f_sp.subs(sol))\n",
    ">\n",
    "(r1, sqrt(-r1^2 + 1)) e^(-1)\n",
    "(r2, -sqrt(-r2^2 + 1)) e^(-1)\n",
    "(0, 0) 0\n",
    "(-1, 0) e^(-1)\n",
    "(1, 0) e^(-1)\n",
    "```\n",
    "Thus the global maxima are along the unit circle (with value $e^{-1}$), and the global minimum is still at $(0,0)$.\n",
    "--->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can also do some calculus involving vectors. Here is the parametrization $r : [0, 2 \\pi] \\to \\mathbb{R}^3$ of a spiral:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "var('t')\n",
    "r=vector([cos(t), sin(t), t]); r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "parametric_plot3d(r, (t, 0, 2*pi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = derivative(r,t); v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The arc length of the path $r$ given by\n",
    "$$\n",
    "\\ell(r) = \\int_{0}^{2 \\pi} | \\dot r(t)| dt.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate(v.norm(),t,0, 2*pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many more interesting functions, and we don't have time to go through all of them. In particular, see [this](https://doc.sagemath.org/html/en/thematic_tutorials/vector_calculus.html) tutorial on vector calculus in SageMath. There, you learn how to define a vector field, compute its divergence and curl, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power series\n",
    "Power series are a very useful part of mathematics. Typically one first encounters them in the form of the Taylor series\n",
    "$$\n",
    "\\sum_{k=0}^\\infty \\frac{f^{(k)}(0)}{k!} t^k\n",
    "$$\n",
    "of a smooth function $f$ around the point $t_0=0$. However, it is also possible to do calculations with more general series\n",
    "$$\n",
    "\\sum_{k=0}^\\infty a_k t^k.\n",
    "$$\n",
    "The easiest way to compute with them in SageMath is to create a ``PowerSeriesRing``, with some formal variable ``t``. Here you have to specify a **base ring**, which is the ring containing the possible coefficients $a_k$ above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R.<t> = PowerSeriesRing(QQ); R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, to compute the Taylor series of the function $f$ at $0$, we can just plug the variable ``t`` of our power series ring into the function ``f``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = exp(t); P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, this computes it to order $t^{19}$, but this can be changed as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R.<t> = PowerSeriesRing(QQ, default_prec=7)\n",
    "P = exp(t); P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that ``default_prec=n`` means the last coefficient that is computed is the one for $t^{n-1}$, which can be somewhat confusing ...\n",
    "\n",
    "The nice thing about power series is that we can add and multiply them, and the Taylor series of $f \\cdot g$ is equal to the product of the series for $f$ and $g$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = sin(t); Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P*Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "taylor(exp(x) * sin(x), x, 0, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the coefficient of $t^k$ of such a power series $P$, we can type ``P[k]``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A cool application of power series are so-called **generating functions**. Given a sequence $(a_k)_{k \\geq 0}$ of real numbers, we say that a function $f(t)$ is a generating function for this sequence, if $\\sum_{k \\geq 0} a_k t^k$ is the Taylor series of $f$ at $t_0 = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "The [Wikipedia page of the Fibonacci numbers](https://en.wikipedia.org/wiki/Fibonacci_number#Power_series) claims that the function\n",
    "$$\n",
    "f(t) = \\frac{t}{1-t-t^2}\n",
    "$$\n",
    "is a generating function for these numbers. Check below if the first coefficients of the Taylor series of $f$ are correct, and use this method to compute the $42$nd Fibonacci number."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution** (uncomment to see)<br>\n",
    "<!---\n",
    "Since we'll later need the coefficient of $t^{42}$, we should better change our power series ring to have a higher precision than $t^{19}$.\n",
    "```\n",
    "R.<t> = PowerSeriesRing(QQ, default_prec=45)\n",
    "f = t/(1-t-t^2); f\n",
    "> t + t^2 + 2*t^3 + 3*t^4 + 5*t^5 + 8*t^6 + 13*t^7 + 21*t^8 + 34*t^9 + 55*t^10 + 89*t^11 + 144*t^12 + 233*t^13 + 377*t^14 + 610*t^15 + 987*t^16 + 1597*t^17 + 2584*t^18 + 4181*t^19 + 6765*t^20 + 10946*t^21 + 17711*t^22 + 28657*t^23 + 46368*t^24 + 75025*t^25 + 121393*t^26 + 196418*t^27 + 317811*t^28 + 514229*t^29 + 832040*t^30 + 1346269*t^31 + 2178309*t^32 + 3524578*t^33 + 5702887*t^34 + 9227465*t^35 + 14930352*t^36 + 24157817*t^37 + 39088169*t^38 + 63245986*t^39 + 102334155*t^40 + 165580141*t^41 + 267914296*t^42 + 433494437*t^43 + 701408733*t^44 + 1134903170*t^45 + O(t^46)\n",
    "```\n",
    "Now we can either read off the coefficients by searching for the ``t^42``-term, or we can extract it as follows:\n",
    "```\n",
    "f[42]\n",
    "> 267914296\n",
    "```\n",
    "Note that a quick check with the pre-implemented ``fibonacci``-function works out:\n",
    "```\n",
    "fibonacci(42)\n",
    "> 267914296\n",
    "```\n",
    "--->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assignments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise (short)\n",
    "Solve Exercise 2 b) from the  [Analysis II-exam](https://www2.math.ethz.ch/education/bachelor/lectures/fs2016/other/analysis2_mavt_matl/pruefungen/pruefung_winter_16.pdf) above.<br>\n",
    "*Remark:* I couldn't get SageMath to solve Exercise 2 a) since it involves some third roots, which are kind of tricky."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution** (uncomment to see)\n",
    "<!---\n",
    "```\n",
    "integrate(1/(x^2-8*x+15),(x,6,Infinity))\n",
    "> 1/2*log(3)\n",
    "```\n",
    "--->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "A partition of a natural number $n \\geq 0$ is one representation of $n$ as an ordered sum of positive integers. For instance, the number $n=5$ has $7$ different partitions:\n",
    "$$\n",
    "5,\\quad 4+1,\\quad 3+2,\\quad 3+1+1,\\quad 2+2+1,\\quad 2+1+1+1,\\quad 1+1+1+1+1\n",
    "$$\n",
    "If $p(n)$ denotes the number of partitions of $n$, then the infinite product\n",
    "$$\n",
    "\\prod_{j=1}^\\infty \\frac{1}{1-t^j} = \\frac{1}{1-t} \\cdot \\frac{1}{1-t^2} \\cdot \\frac{1}{1-t^3} \\cdots\n",
    "$$\n",
    "is a generating function of $p(n)$. \n",
    "* (theoretical) Use the formula for the geometric series to show that the finite product\n",
    "$$\n",
    "\\prod_{j=1}^m \\frac{1}{1-t^j} = \\frac{1}{1-t} \\cdots \\frac{1}{1-t^m}\n",
    "$$\n",
    "is a generating series for the number $p(n,m)$ of partitions of $n$ with maximal part at most $m$.<br>\n",
    "In particular, since for $m \\geq n$ we have $p(n,m)=p(n)$, this explains why the infinite product above is a generating series for $p(n)$.\n",
    "* Use this knowledge to write a function ``partition_number(n, m=None)`` to compute the numbers $p(n,m)$.<br>\n",
    "*Remark:* In the above case, it is useful to set the optional parameter ``m`` to the Python value ``None`` as default. Then, at the beginning of the function, you should check whether ``m`` is ``None`` and if yes give it a reasonable value.\n",
    "* Check your function in a few cases (e.g. $n=5$ using the enumeration of partitions of $5$ above)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solution** (uncomment to see)<br>\n",
    "<!---\n",
    "* For an explanation of the first statement, see e.g. the following [Wikipedia page](https://en.wikipedia.org/wiki/Partition_function_(number_theory)#Generating_function).\n",
    "* Here is a possible implementation. Note that we have to set ``default_prec=n+1`` to be able to extract the coefficient of $t^n$ later.\n",
    "\n",
    "```\n",
    "def partition_number(n, m=None):\n",
    "    if m is None:\n",
    "        m = n\n",
    "    R.<t> = PowerSeriesRing(QQ,default_prec=n+1)\n",
    "    f = prod([1/(1-t^j) for j in range(1,m+1)])\n",
    "    return f[n]\n",
    "```\n",
    "* We check a few cases:\n",
    "```\n",
    "partition_number(5)\n",
    "> 7\n",
    "partition_number(5,3)\n",
    "> 5\n",
    "```\n",
    "--->"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.1",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}